
****************************;
* Pap model - table 2 PAP;
****************************;
proc import datafile = 'yourpathway\pap_wqualclin.csv'
	out = pap_1 
	dbms = csv;
	guessingrows = max;
	getnames = yes;
run;
data analysis;
set pap_1;
month=month(month_start);
year=year(month_start);
if year=2020 then month+12;
if year=2021 then month+24;
if month<15 then do;
period="1_Pre Pandemic";
	month_int=month;
	covid=0;
end;
if month>=15 then do;
period="2_COVID";
	month_int=month;
	covid=1;
end;
if month = 15 then pap_patient_tests = .;
run;
proc glimmix data=analysis;
class delivery_site_id_masked covid(ref=first);
model pap_patient_tests= covid month_int covid*month_int /
link=log offset = offset_denom d=poisson solution CL;
random intercept /subject=delivery_site_id_masked;
ods output parameterestimates=parameters;
/*out predicted values in case you want them*/
output out=predicted pred(blup ilink)=predicted
	lcl(blup ilink)=lower
	ucl(blup ilink)=upper;
run;
*exp in paper;
data parameters1; set parameters;
est1 = exp(Estimate);
l1 = exp(lower);
u1 = exp(upper);
run;
***********************************;
* PAP RATE - figure 2 panel C;
ODS OUTPUT lsmeans=lsm1(KEEP = month_start estimate lower upper);
proc genmod data=pap_1;
	class delivery_site_id_masked month_start(ref = "2019-01-01") ;
	model pap_patient_tests = month_start
	/ dist=poisson link=log offset=offset_denom type3;
	repeated subject=delivery_site_id_masked / type=ar(1);
	lsmeans month_start  / cl;
run;
data Lsm1; set Lsm1;
estimate1 = exp(estimate);
lower1 = exp(lower);
upper1 = exp(upper);
run;
proc sort data = lsm1;
by month_start;
run;
ods layout gridded columns=1;
*ods region;
ods graphics on / width=7in;
ods graphics on / height=4in;
 ods html image_dpi=400;
proc sgplot data=lsm1 noautolegend;     
styleattrs datacontrastcolors=(gold navy); 
   scatter x=month_start y=estimate1 / 
						   yerrorlower=lower1                                                                                            
                           yerrorupper=upper1   
						   errorbarattrs =(color=black)
                           markerattrs=(color = black symbol=CircleFilled)
						  ;   
 
   series x=month_start y=estimate1 /  lineattrs=( color = black pattern=1);  
   refline '2020-03-01' / axis= x lineattrs=(thickness=3 color=black pattern=dash);
   *reg  x=month_start y=estimate1 /  lineattrs=( color = black pattern=1) markerattrs = (color = red);
   *keylegend  "eighty"/ location = inside 
   						position = bottomright;
   xaxis label = "Date"  interval = month  display=(nolabel);
   yaxis label = "Rate of PAP Tests per 1000 Eligible Patients"; 
run;  
***********************************;
*plot up to date Figure 1 panel B;
proc sort data = pap_1;
by month_start;
run;
proc means data = pap_1;
by month_start;
var pct_utd_pap_uds1;
output out = total MEAN=mean
STDERR=stderr LCLM=lclm Uclm=uclm ;
run;
data total1; set total;
length delivery_site_id_masked $12;
delivery_site_id_masked = "Overall";
run;
proc sort data = total1;
by delivery_site_id_masked ;
run;
data pap_spec;
set pap_1;
	if delivery_site_id_masked in ("2030_02", "6880_02", "2053_01", "6008_01", "7778_02", "3338_01",
           "4222_01", "8008_02", "1118_01", "5442_01", "7201_03",  "3687_01") 

	then output;
length Sample $10.;
Sample = "Recruited";
length pct_utd_pap_uds1 8.;
pct_utd_pap_uds1 = input(pct_utd_pap_uds, 8.);
run;
proc sort data = pap_spec;
by delivery_site_id_masked ;
run;
data together;
merge pap_spec total1;
by delivery_site_id_masked ;
run;
data together1;
set together;
if delivery_site_id_masked = "Overall"
       then overall_total = mean;
else clinic_total = pct_utd_pap_uds1;
label overall_total="Overall Average" clinic_total="Clinics" ;
run;
proc sort data=together1;
by delivery_site_id_masked;
run;
data together1;
set together1;
by delivery_site_id_masked;
retain label 0;
if first.delivery_site_id_masked then label=label+1;
if label = 6 or label = 2 or label = 3 or label = 7 or label = 8  then clinic1_total = clinic_total;
else clinic1_total = .;
if label = 6 or label = 2 or label = 3 or label = 7 or label = 8 then clinic2_total = .;
else clinic2_total = clinic_total;
run;

ods graphics on / width=7in;
ods graphics on / height=4in;
ods html image_dpi=400;
proc sgplot data = together1 noautolegend nowall noborder;
styleattrs datacontrastcolors=(mediumblue orangered gold purple skyblue brown powderblue coral Green fuchsia red thistle);

	series x = month_start y = clinic2_total / 
											group = label CURVELABEL curvelabelloc=outside 
											CURVELABELATTRS=(size=9.5) lineattrs=(pattern=solid);
	series x = month_start y = clinic1_total/ group = label CURVELABEL curvelabelloc=inside CURVELABELPOS = start 
											CURVELABELATTRS=(size=8.5) lineattrs=(pattern=solid);
	series x = month_start y = overall_total / group = delivery_site_id_masked lineattrs=(pattern=dash color = black thickness = 4);
	xaxis label = "Date"  interval = month display=(nolabel);
	refline 2020-03-01 / axis = x lineattrs=(color = black pattern = solid);
    yaxis min = 0 max = 100 label = "Percent Up to Date for Cervical Cancer Screenings"; 
run;
*****************************************;
**** CRC/Colorectal Cancer Screening ****;
*****************************************
flexsig_patient_tests 
fobt_patient_tests 
fit_patient_tests 
colonoscopy_patient_tests;
proc import datafile = 'yourpathway\colo.csv'
	out = colo 
	dbms = csv;
	guessingrows = max;
	getnames = yes;
run;
data colo_1;
set colo;
keep delivery_site_id_masked month_start pap_patient_tests 
offset_denom colonoscopy_patient_tests combo_f pct_utd_crc_uds1;
run;

data analysis_colo;
set colo_1;
month=month(month_start);
year=year(month_start);
if year=2020 then month+12;
if year=2021 then month+24;
if month<15 then do;
period="1_Pre Pandemic";
	month_int=month;
	covid=0;
end;
if month>=15 then do;
period="2_COVID";
	month_int=month;
	covid=1;
end;
if month = 15 then pap_patient_tests = .;
run;
*table 2 colo;
proc glimmix data=analysis_colo;
class delivery_site_id_masked covid(ref=first);
model colonoscopy_patient_tests= covid month_int covid*month_int /
link=log offset = offset_denom d=poisson solution CL;
random intercept /subject=delivery_site_id_masked;
ods output parameterestimates=parameters;
/*out predicted values in case you want them*/
output out=predicted pred(blup ilink)=predicted
	lcl(blup ilink)=lower
	ucl(blup ilink)=upper;
run;
data parameters; set parameters;
est1 = exp(Estimate);
l1 = exp(lower);
u1 = exp(upper);
run;
*table 2 fit/fobt;
proc glimmix data=analysis_colo;
class delivery_site_id_masked covid(ref=first);
model combo_f= covid month_int covid*month_int /
link=log offset = offset_denom d=poisson solution CL;
random intercept /subject=delivery_site_id_masked;
ods output parameterestimates=parameters_fit;
/*out predicted values in case you want them*/
output out=predicted pred(blup ilink)=predicted
	lcl(blup ilink)=lower
	ucl(blup ilink)=upper;
run;
data parameters_fit; set parameters_fit;
est1 = exp(Estimate);
l1 = exp(lower);
u1 = exp(upper);
run;
ODS OUTPUT lsmeans=lsm1(KEEP = month_start1 estimate lower upper);
proc genmod data=colo;
	class delivery_site_id_masked month_start1(ref = "2019-01-01") ;
	model colonoscopy_patient_tests = month_start1
	/ dist=poisson link=log offset=offset_denom type3;
	repeated subject=delivery_site_id_masked / type=ar(1);
	lsmeans month_start1  / cl;
run;

ODS OUTPUT lsmeans=lsm4(KEEP = month_start1 estimate lower upper);
proc genmod data=colo;
	class delivery_site_id_masked month_start1(ref = "2019-01-01") ;
	model combo_f = month_start1
	/ dist=poisson link=log offset=offset_denom type3;
	repeated subject=delivery_site_id_masked / type=ar(1);
	lsmeans month_start1  / cl;
run;
data Lsm1; set Lsm1;
estimate1 = exp(estimate)*1000;
lower1 = exp(lower)*1000;
upper1 = exp(upper)*1000;
run;
data Lsm4; set Lsm4;
estimate1 = exp(estimate)*1000;
lower1 = exp(lower)*1000;
upper1 = exp(upper)*1000;
run;
proc sort data = lsm1;
by month_start1 ;
run;
proc sort data = lsm4;
by month_start1 ;
run;
*figure 2 panel A - colo rate;
proc sgplot data=Lsm1 noautolegend ;     
styleattrs datacontrastcolors=(gold ); 
   scatter x=month_start1 y=estimate1 / 
						   yerrorlower=lower1                                                                                            
                           yerrorupper=upper1   
						   errorbarattrs =(color=CX03744B)
                           markerattrs=(color = CX03744B symbol=CircleFilled)
						  ;   
 
   series x=month_start1 y=estimate1 /  lineattrs=( color = CX03744B pattern=1);
	  * reg  x=month_start y=estimate1 /  lineattrs=( color = purple pattern=1);
 
   *keylegend  "eighty"/ location = inside 
   						position = bottomright;
   xaxis label = "Date"  interval = month  display=(nolabel);
   yaxis label = "Rate of Colonoscopies per 1000 Eligible Patients"; 
run;
*********************************
*figure 2 panel b -fobt/fit rate;
proc sgplot data=lsm4 noautolegend;     
*styleattrs datacontrastcolors=(gold ); 
   scatter x=month_start1 y=estimate1 / 
						   yerrorlower=lower1                                                                                            
                           yerrorupper=upper1   
						   errorbarattrs =(color=navy)
                           markerattrs=(color = navy symbol=CircleFilled)
						  ;   
 
   series x=month_start1 y=estimate1 /  lineattrs=( color = navy pattern=1);  
   *reg  x=month_start y=estimate1 /  lineattrs=( color = purple pattern=1);
 
   *keylegend  "eighty"/ location = inside 
   						position = bottomright;
   xaxis label = "Date" interval = month  display=(nolabel);
   yaxis label = "Rate of FOBT/FIT Tests per 1000 Eligible Patients"; 
run; 
***************************************************;
************Colo Percent up to date****************;
***************************************************;
* colo up to date ;
proc sort data = colo_1;
by month_start;
run;

proc means data = colo_1 ;
by month_start;
var pct_utd_crc_uds1;
output out = total MEAN=mean
STDERR=stderr LCLM=lclm Uclm=uclm ;
run;
data total1; set total;
length delivery_site_id_masked $12;
delivery_site_id_masked = "Overall";
*if _STAT_ = "MEAN" and month_start ^= . then output;
run;
proc sort data = total1;
by delivery_site_id_masked ;
run;
data colo_spec;
set colo_1;
length Sample $10.;
Sample = "Recruited";
	if delivery_site_id_masked in ("2030_02", "6880_02", "2053_01", "6008_01", "7778_02", "3338_01",
           "4222_01", "8008_02", "1118_01", "5442_01", "7201_03",  "3687_01") 
	then output;
run;
proc sort data = colo_spec;
by delivery_site_id_masked ;
run;
data together;
merge colo_spec total1;
by delivery_site_id_masked ;
run;

*plot up to date;
data together1;
set together;
if delivery_site_id_masked = "Overall"
       then overall_total = mean;
else clinic_total = pct_utd_crc_uds1;
label overall_total="Overall Average" clinic_total="Clinics" ;
run;
proc sort data=together1;
by delivery_site_id_masked;
run;
data together1;
set together1;
by delivery_site_id_masked;
retain label 0;
if first.delivery_site_id_masked then label=label+1;
*if label = 6 or label = 2 or label = 3 or label = 7 or label = 8  then clinic1_total = clinic_total;
*else clinic1_total = .;
*if label = 6 or label = 2 or label = 3 or label = 7 or label = 8 then clinic2_total = .;
*else clinic2_total = clinic_total;
run;

 ods html image_dpi=400;
proc sgplot data = together1  nowall noborder;
	series x = month_start y = clinic_total / group = label CURVELABEL curvelabelloc=outside 
											  CURVELABELATTRS=(size=9.5) lineattrs=(pattern=solid);
	series x = month_start y = clinic_total / group = delivery_site_id_masked ;
	series x = month_start y = overall_total / group = delivery_site_id_masked lineattrs=(pattern=dash color = black thickness = 4);
	xaxis label = "Date" interval = month  display=(nolabel);
    yaxis min = 0 max = 100 label = "Percent Up to Date for Colorectal Cancer Screenings"; 

run;
